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Abstract. In this paper, we propose a discretization for the (nonlinearized) 
compressible Stokes problem with a linear equation of state p = p, based on 
Crouzeix-Raviart elements. The approximation of the momentum balance is 
obtained by usual finite element techniques. Since the pressure is piecewise 
constant, the discrete mass balance takes the form of a finite volume scheme, in 
which we introduce an upwinding of the density, together with two additional 
stabilization terms. We prove a priori estimates for the discrete solution, which 
yields its existence by a topological degree argument, and then the convergence 
of the scheme to a solution of the continuous problem. 



1. Introduction 

The problem addressed in this paper is the system of the so-called barotropic 
compressible Stokes equations, which reads: 

(1.1a) - Au + Vp = f 

(1.1b) div(p«)=0 

(l.lc) p = g(p) 

where p, u and p stand for the density, velocity and pressure in the flow, respec- 
tively, and / is a forcing term. The function g is the equation of state used for 
the modelling of the particular flow at hand, which may be the actual equation of 
state of the fluid or may result from assumptions concerning the flow. Here, we 
only consider the following relation, which corresponds to an isothermal flow of a 
perfect gas: 

(1.2) e(p)-A P 

where A is a positive constant. Since the sound velocity is defined by c 2 = dp/dp, 
A = Msl 2 /V 2 , where Ma is the Mach number and V is a characteristic velocity of 
the flow. This system of equations is posed over fi, a bounded domain of M. d , d < 3 
supposed to be polygonal (d = 2) or polyhedral (d — 3). It is supplemented by 
homogeneous boundary conditions for u, and by prescribing the total mass M of 
the fluid: 

(1.3) / pdx = M 

Jn 
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where M is a positive real number. 

In this paper, we study a numerical scheme for the solution of this problem, 
which combines low order finite element and finite volume techniques, and is very 
close to a scheme which was proposed for the solution of barotropic Navier-Stokes 
equations in |10j and further extended to two-phase flows in the resulting 

code is today currently used at the French Institut de Radioprotection et de Surete 
Nucleaire (IRSN) for "real-life" studies in the nuclear safety field. Up to now, 
stability (in the sense of conservation of the entropy) is known for these schemes, 
and numerical experiments show convergence rates close to one in natural energy 
norms. Our goal is now to prove their convergence. This work is the first one in 
this direction, and we address here the probably simplest toy problem, restricting 
ourselves to the steady case, to creeping flows (i.e. omitting the convection term in 
the momentum balance equation) and to a linear equation of state. The extension 
to laws where p varies linearly with p 1 ^ 1 , where 7 > 1 is a coefficient which is 
specific to the considered fluid, which are typically obtained for isentropic flows 
of perfect gases, is the object of a further paper (part II of the present one); the 
additional difficulty posed by this further study is to prove the strong convergence 
for the density, which necessitates to adapt P.L. Lions' "effective pressure trick" 
[16j at the discrete level. Finally, for the sake of simplicity, we use here a simplified 
form of the diffusion term (-Am) but it is clear from the subsequent developements 
that the presented theory holds for any linear elliptic operator (and in particular 
for the usual form of the viscous term for compressible constant viscosities flows). 

The finite element - finite volume discretization which is chosen here is moti- 
vated by the fact that we wish the approximate density to be positive, as in the 
continuous model, in order to be compatible with the physics. Moreover, the proof 
of convergence of a numerical approximation to (jl.ip requires estimates on both 
velocity and pressure or density, and the density positivity is very useful to obtain 
these estimates. A classical way to ensure positivity is to use a finite volume up- 
winding technique in the discretization of the term div(pu). This technique is easily 
set up if the discrete velocities are located on the edges and densities and pressures 
at the cell centres, which is the reason why we choose the Crouzeix-Raviart finite 
elements for the velocities and cell centred finite volumes for the densities. 

This paper is organized as follows. The discretization is first described (section 
[2]) , and we prove an L 2 compactness result for sequences of Crouzeix-Raviart func- 
tions with bounded broken H 1 semi-norm (section [3]) . Then the proposed scheme 
is given (section |4]), and the above-mentioned compactness result yields the conver- 
gence of (sub-)sequences of discrete solutions to a limit, thanks to a priori estimates 
which are given in section [5] Finally, this limit is shown to be a solution to the 
continuous problem in section [6l 

To our knowledge, this convergence proof is the first one for the genuine (non- 
linear) compressible Stokes problem; a linearized version of this system is adressed 
in previous works [H EH HI HH □ • 

2. Discrete spaces and relevant lemmata 

Let T be a decomposition of the domain in simplices. By £(K), we denote 
the set of the edges (d = 2) or faces (d = 3) a of the element K e T; for short, 
each edge or face will be called an edge hereafter. The set of all edges of the mesh 
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is denoted by £ ; the set of edges included in the boundary of O is denoted by £ ey ± 
and the set of internal ones {i.e. £ \ £ ex t) is denoted by £ ln t- The decomposition T 
is supposed to be regular in the usual sense of the finite element literature {e.g. [3]), 
and, in particular, T satisfies the following properties: Cl = {_} Ke q- K; if K, LeT, 
then A > nL = 0,A v nL = 0isa vertex of the mesh or K n L is a common edge of K 
and L, which is denoted by K\L. For each internal edge of the mesh a = K\L, tikl 
stands for the normal vector of cr, oriented from K to L (so that tikl = -mlk). By 
\K\ and \a\ we denote the measure, respectively, of the element K and of the edge 
<7, and Kk and h a stand for the diameter of K and a, respectively. We measure 
the regularity of the mesh through the parameter 6 defined by: 

(2.1) 6 = inf K eT}U{^,^; * = K\L g & mt }, 

where stands for the diameter of the largest ball included in K. Note that the 
following inequality holds: 

(2.2) h a \a\ < 2 0- d \K\, VK g T, Vct g S(A'). 

Indeed, this relation is derived by noting that | cr] < h d K and A > with 
c = 7r/4 in 2D and c = 7r/6 in 3D; it will be used throughout this paper. Finally, 
as usual, we denote by h the quantity maxxgr Hk- 

The space discretization relies on the Crouzeix-Raviart element (see [4] for the 
seminal paper and, for instance, p. 199-201] for a synthetic presentation). The 
reference element is the unit <i-simplex and the discrete functional space is the space 
Pi of affine polynomials. The degrees of freedom are determined by the following 
set of nodal functionals: 

(2.3) {F a , erg F a (v) = \tr\~ 1 f «d 7 . 

The mapping from the reference element to the actual one is the standard affine 
mapping. Finally, the continuity of the average value of the discrete functions {i.e., 
for any function v, F a {v)) across each face of the mesh is required, thus the discrete 
space Vh is defined as follows: 

V h = {»£ L 2 {VL) : v\ K g Pi (if), VA g T; F a {v) continuous 

(2.4) 

across each edge a G £i nt ; F a {v) = 0, Vcr 6 £ e xt}- 

The space of approximation for the velocity is the space Wh of vector valued func- 
tions each component of which belongs to Vh- Wh = {Vh) d . The pressure is ap- 
proximated by the space Lh of piecewise constant functions: 

L h = {q e L 2 {Q) : q\ K = constant, VK E T} . 

Since only the continuity of the integral over each edge of the mesh is imposed, 
the functions of Vh are discontinuous through each edge; the discretization is thus 
nonconforming in ii 1 (£!) d . We then define, for 1 < i < d and v g Vh, dh,i v as the 
function of L 2 (f2) which is equal to the (piecewise constant) derivative of v with 
respect to the i th space variable almost everywhere. This notation allows to define 
the discrete gradient, denoted by V/j, for both scalar and vector valued discrete 
functions and the discrete divergence of vector valued discrete functions, denoted 
by div h . 
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The Crouzeix-Raviart pair of approximation spaces for the velocity and the pres- 
sure is inf-sup stable, in the usual sense for "piecewise H 1 " discrete velocities, i.e. 
there exists c; > independent of the mesh such that: 



J K Q dWvdx [ qdxv h vdx 

_ . J o 



VqeL h , sup j, = sup > Q \\q - g m || L =(n) 

vew h \\ v \h,b vew h \\ v h,b 

where g m is the mean value of q over Q and || • Hi^ stands for the broken Sobolev 
H 1 semi-norm, which is defined for any function v £ Vh or v S Wh by: 

ll"lli,6 = E / |V«| 2 dx= / |V^| 2 dx. 



a 



This broken Sobolev semi-norm is known to control the L 2 norm by an extended 
Poincare inequality 19, proposition 4.13], in the sense that for any function v G Vh, 
IM|i,6 5= c ll u l|L 2 (n) where the real number c only depends on the computational 
domain. 

We also define a discrete semi- norm on Lh, similar to the usual finite volume 
discrete H 1 semi-norm, weighted by a mesh-dependent coefficient: 

V9GL,, \q\l = {h K + h L f M (q K -q L ) 2 . 

a=K\L 

From the definition (|2.3[) , each velocity degree of freedom can be associated to 
an element edge. Consequently, the velocity degrees of freedom are indexed by the 
number of the component and the associated edge, thus the set of velocity degrees 
of freedom reads: 

a G £ int , 1 < i < d}. 
We denote by (f> a the usual Crouzeix-Raviart shape function associated to er, i.e. 
the scalar function of Vh such that F a (4> a ) = 1 and F<,>((f) a ) — 0, W G £ \ {a}. 

Similarly, each degree of freedom for the pressure is associated to a cell K, and 
the set of pressure degrees of freedom is denoted by {pk, K G T}. 

Finally, we define by the following interpolation operator: 



(2.5) 



Hj(fi) — » V h 



' i ^ r h v = 

This operator naturally extends to vector- valued functions (i.e. to perform the 
interpolation from Hj(f2) d to Wh), and we keep the same notation for both the 
scalar and vector case. The properties of r/j are gathered in the following lemma. 
They are proven in [4]. 

Lemma 2.1. Let 8q > and let T be a triangulation of the computational domain 
Q such that 9 > 9q, where 9 is defined by (|2.1jl . The interpolation operator rh 
enjoys the following properties: 

(1) Preservation of the divergence: 

(2.6) Vv G Hl(fl) d , Vq G Lh, / q div h (r h v)dx = / qdivvdx. 

Jn Jn 
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(2) Stability: 

(2.7) WeHS(n), \\r h v\\ ltb <ci(0 o )M H i(n)- 

(3) Approximation properties: 

Vv g H 2 (fflnHj(0), vjf g r, 

(2.8) MM 2 

F - r h v\\ L 2 {K) +h K \\V h (v-r h v)\\ L 2 (K) < c 2 (9 ) h K \v\ H * {K) . 

In both above inequalities, the notation Ci(#o) means that the real number a 
only depends on 9$, and, in particular, not on the parameter h characterizing the 
size of the cells; this notation will be kept throughout the paper. 

The following lemma is known {e.g. [5J lemma 3.32]); we give its (elementary) 
proof for the sake of completeness. 

Lemma 2.2. Let 9q > and let T be a triangulation of the computational domain 
Q such that 9 > 9q, where 9 is defined by (|2.ip . and Vh be the space of Crouzeix- 
Raviart discrete functions associated to T, as defined by (|2.4p . Then there exists a 
real number c(9q) such that the following bound holds for any v £ Vh-' 



Er / M 2 dj<c(e )\\v\\l b , 



<re£ 

where, on any a G & m t, M stands for the jump of v across a and, on any a G £ ex t, 
[v] = v. 

Proof. For any control volume K of the mesh, we denote by (Vi?)k the (constant) 
gradient of the restriction of v to K . With this notation, using the continuity of v 
across a at the mass center x a of any internal edge and the fact that v vanishes at 
the mass center x a of any external edge, we get: 

£_L f[ v fd 7 = Y f {{{Vv) K -{Vv) L )-{ X -x a )f ^ 

£ ^J^v) K -(x-x a )) 2 d 7 . 



•yes aeSi 

o=K\L 



We thus have: 

Er / H 2 d 7 <2 £ (\(Vv) K \ 2 + \(Vv) L f)+ J2 K\a\\{Vv) K \\ 

a=K\L <j£E{K) 

and the result follows by regularity of the mesh. □ 
The proof of the following trace lemma can be found in [211 section 3] . 

Lemma 2.3. Let T be a given triangulation offl and K be a control volume ofT, 
h,K its diameter and a one of its edges. Let v be a function o/H 1 (if). Then the 
following inequality holds: 

IM|l=(<7) < [d J ( H|l=(K) + h K ||Vu|| L 2 (K) ) . 
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We will also need the following Poincare ineqality: 



(2.9) 



g T, Vv e H\K), \\v - v m , K \\ h 2 (K) < - h K ||Vw|| L 2 (K) 



where stands for the mean value of v over K . This relation is proven for any 
convex domain in 181. 



We are now in position to prove the following technical result. 

Lemma 2.4. Let 8q > and let T be a triangulation of the computational domain 
f2 such that > 6q, where 9 is defined by (|2.1[) ; let {a a ) a ^£ iTit be a family of real 
numbers such that Vcr G £ m t > a a < 1 and let v be a function of the Crouzeix-Raviart 
space Vh associated to T . Then the following bound holds: 



E 

crGfint 



a a H/d7 



<c(6 )h\\v\\ 1 . b |/| H i(n) , V/ G Hj(n). 



where the real number c(0q) only depends on 8q and on the domain fl. 

Proof. Since the integral of the jump across any edge of the mesh of a function of 
Vh is zero, we have, for any a G £i n t- 



a<r [v] /d7 = / a a [v] (f - f a ) fry, 



where f a is any real number. Using the Cauchy-Schwarz inequality, first in L 2 (ct) 
then in M card ( £ ) we thus get: 



E 



a CT [v] f dj 



< E 



H 2 d 7 



1/2 



(/-/ CT ) 2 d7 



1/2 



< 



i 1/2 r 

E T- /H 2 d 7 E hj(f-f a ffr f 



1/2 



By lemma l2~2| the first term of the latter product is bounded by c(0q) • For 

the second one, choosing arbitrarily one adjacent simplex to each edge and applying 
the above trace lemma [231 we get: 



n < ]r 2dK M ( ||/ _ fa \\i 2{K) + hi ||v/ii 2 2 



(K) 



(<ree(K)) 



Choosing for / CT the mean value of / on K and using (|2.9|) . we thus get: 
2?< 2 2d(l + ^)^M^|| V /|| 2 s 



fGfint 

(aG£(if)) 



lL 2 (if) 



and the result follows by observing that the H 1 semi-norm of / on K appears at 
most (d + 1) times in the summation and using the regularity of the mesh. □ 
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3. A COMPACTNESS RESULT 

The aim of this section is to state and prove a compactness result for || ■ Hi^ 
bounded sequences of discrete functions. We begin by a preliminary lemma. 

Lemma 3.1. Let 9q > and let T be a triangulation of the computational domain 
such that 9 > 9q, where 9 is defined by (|2.ip ; for a G £, let Xa be the function 
defined by: 

Xa ■ K d x R d — ► {0, 1} 

(x, y) i ► Xa(x, y) = 1 if [x, y] n a ^ 0, x„{x, y) = otherwise, 

where x and y are two points of M. d . Then there exists a family of positive real 
numbers (d a )a-eE such that: 

(1) for any a G £, d a = C\{6q) h a , 

(2) for any points x and y of M. d (possibly located outside Vt), the following 
inequality holds: 

^2xa{x 7 y) d a < c 2 (6» ) (\y-x\ + h) 

Proof. We first deal with the two-dimensional case and with quasi-uniform meshes 
(i.e. the bound we first prove blows up when maxK£T{h/hK) tends to infinity). 

Let T be a triangulation of a two-dimensional domain J7, K a triangular cell 
of T and a an edge of K. Without loss of generality, we suppose that a is the 
segment (0, h a ) x and we denote by £k the diameter of the largest ball included 
in K and by hx the diameter of K. We denote by z a the opposite vertex to er; 
the first coordinate of z a is necessarily lower than hpc while its second coordinate 
is necessarily greater than (in the opposite case, no ball of diameter £k would 
be included in K). It thus follows (see figure [T|): 

(1) that the rectangular domain uj a = (h a /3,2h a /3) x (0, h^K is 
included in K, 

(2) that, if the similar construction is performed for another edge a' of K to 
obtain av, uj a and lo' g do not intersect. 

We denote by d a the quantity d a = h a £,K I '(12/&K-) • We thus have d a > (9/12) h a , 
where 9 is the parameter defined by 1 2 . H 

We now perform this construction for each edge a of the mesh. If a G £ xt, there 
is only one possible choice for K (the adjacent cell to a); if cr G £i n t, cr = K\L, we 
choose either K or L. Let x and y be two points of M 2 . Let t( x ,y) be the vector 
given by: 

y — x 
t{x ^ = ]y—^\ 

and ni x ,y) a normal vector to t(~ ty y We denote by S^ x y ^ the rectangle defined by: 

S( X , V ) = {x + £i t( x>y ) + 6 ri( Xt y), 6 e (-ft, |y - x\ + h), ^ 2 g (-h, +h)} 

For each edge intersected by the segment [x,y] (i.e. for each edge a such that 
X<r{x,y) = 1), the rectangle uj a is included in Sr^y)] thus, since these domains Lo a 
and are disjoint: 
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Figure 1 . Notations for the control volume K 



and thus: 

X'fo y)\d a K<2h (\y -x\ + 2h), 

which concludes the proof if max#- 6 7-(/i//ix) is supposed to be bounded. 

The extension to the three-dimensional case is straightforward, since it only 
necessitates to adapt the construction of the domains Lo a . Finally, giving up the 
assumption that maxKeTih/hx) is bounded only needs a more careful definition 
of the domain S^ x y ^, replacing the parameter h by a local value. □ 



The following bound provides an estimate of the translates of a discrete function 
v as a function of . 

Lemma 3.2. Let 8q > and let T be a triangulation of the computational domain 
f2 such that 9 > 6q, where 8 is defined by (|2.1[) ; let be the space of Crouzeix- 
Raviart discrete functions associated to T, as defined by (|2.4[) . Let v be a function 
ofVh; we denote by v the extension by zero of v to K d . Then the following estimate 
holds: 

Vr, e R d , IKK" +V)~ 5(-)IIl=.(r*) < c(ffo) \V\ (\v\ + h) \\v\\l b . 



Proof. We follow the proof of a similar result for piecewise constant functions, 
namely Lemma 9.3, pp. 770-772]. Let r] G K d be given, v be a Crouzeix-Raviart 
discrete function and v its extension by zero to K d . With the definition of the 
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function Xa of Lemma [37X1 the following identity holds for any x € M d : 

v(x + rj) - v(x) = ^2xa(x,x + r]) [v\{y x _ VM ) + / V h v(x + srj) ■ r/ds 



T x {x) 



T 2 (x) 



where y x v a stands for the intersection between the line issued from x and of 
direction r\ and the hyperplane containing a. Defining for each edge a of the mesh 
a real positive number d a such that Lemma 13.11 holds, by the Cauchy-Schwarz 
inequality, we get for T\(x): 

{Ti{x)f < fexa(x,x + »|) ^^'^ (jZxAx.x + rj) A 

Integrating now over Mr, we thus obtain: 

/ mx)) 2 dx < c 2 (9 ) (\v\ + h)[Jl[ Xa(x,x + V ) M(y . w)2 dx] 

Let Qa,ri = {x = y + srf.y 6 a and s € [—1,0]}. Noting that the function x 
Xa(x,x + rj) is in fact the characteristic function of Q a .ri, we get that: 

/ Xa(x,x + r]) {[v]{y x ^^)) 2 dx = I (M(y m ^)) 2 dx 

= \n a -v\J J {[v]{y)f dy ds, 
where n CT is a unit normal vector to a. Therefore: 

/ (T^xtfdxKcM (\r,\ + h) \ti\J2t [ (H(y)) 2 dy, 

and thus, by choice of d a : 



(3.1) (Ti(x)r 

f'l NTn 

On the other hand, by the Cauchy-Schwarz inequality, we have for 



(^(a;)) 2 < |77| 2 f \V h v(x + s V )\ 2 ds, 
Jo 



and thus, using the Fubini theorem and remarking that V/jW vanishes outside Q, 
we get: 

(3.2) / (T 2 (x)) 2 dx<H 2 ||< b . 

The result then follows thanks to the inequality \v(x + 77) — v(x)\ 2 < 2(Ti(x)) 2 + 
2(T 2 (x)) 2 , to the bounds (EH) and $£2$ and to Lemma [2^1 □ 

We are now in position to state the following compactness result. 

Theorem 3.3. Let («' m ') m eN be a sequence of functions satisfying the following 
assumptions: 
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(1) Vra E N, there exists a triangulation of the domain T^" 1 ' such that t>( m ) £ 
V^ m \ where V^" 1 ' 1 is the Crouzeix-Raviart space associated to 7~( m ) defined 
by (12. 4p . and the parameter #( m ) characterizing the regularity of T^ m > is 
bounded away from zero independently of m, 

(2) the sequence (v^ m ') me ^ * s uniformly bounded with respect to the broken 
Sobolev H 1 semi-norm, i.e.; 

Vra G N, \\v (m) \\i,b < C 
where the real number C does not depend on m and \\ ■ stands for the 
broken Sobolev H 1 semi-norm associated to 7~( m ) (with a slight abuse of 
notation, namely dropping for short the index ty7n ^ pointing the dependence 
of the norm with respect to the mesh). 

Then, possibly up to the extraction of a subsequence, the sequence («' m ') m gB con- 
verges strongly in L 2 (f2) to a limit v such that v G Hj(f2). 



Proof. The result follows from the translates estimates of lemma 13.21 The com- 
pactness in L 2 (f}) of the sequence is a consequence of the Kolmogorov theorem (see 
e.g. [6] theorem 14.1, p. 833] for a statement of this result). The fact that the 
limit belongs to Hj(J7) follows from the particular expression for the bound of the 
translates and is proven in [SJ theorem 14.2, pp. 833-834]. □ 

4. The numerical scheme 

Let p* be the mean density, i.e. p* = M/|f2| where |f2| stands for the measure of 
the domain fl. We consider the following numerical scheme for the discretization 
of Problem ([L~T]) : 

(4.1a) Vu G Wh, / Vft,w : V^u da; — / pdivhvdx = / f vdx, 
Jn Jn Jn 

(4.1b) VKeT, (<k QM - v"k QiPL)) + h a \K\ (g(p K ) - p*) 

a=K\L „ 

J stab, 1 

+ {hK + h^Y \q{pk) + q{pl)\ (q(pk) - q{j>l)) = 0, 

a=K\L " 



?stab,2 

where v+ K and v~ K stands respectively for v~^ K = max(v CT x, 0) and v~ K = 
— mm(v a ^x, 0) with v a _K = \o~\u a ■ tlkl = v ai< ~ v a i<- Note that the up- 
winded convection term Y^o-=K\L ( v <t k Q(Pk) — v ^k q(Pl)j may a l so be written: 

Y,cr=K\L V ^,KPa, with 

(4.2) p fj = 



PK if V<r,K > 0, 

Pl otherwise. 



Equation (|4.1a[) may be considered as the standard finite element discretization of 
(|l.lap . Since the pressure is piecewise constant, the finite element discretization of 
(|l.lb[) . i.e. the mass balance, is similar to a finite volume formulation, in which we 
introduce the standard first-order upwinding and two stabilizing terms. The first 
one, i.e. T sta b lj guarantees that the integral of the density over the computational 



A FE-FV SCHEME FOR THE ISOTHERMAL COMPRESSIBLE STOKES PROBLEM 11 



domain is always M (this can easily be seen by summing the second relation for K G 
T). The second one, i.e. T s t a b,2, is useful in the convergence analysis. It may be seen 
as a finite volume analogue of a continuous term of the form div (|p| Vp) weighted by 
a mesh-dependent coefficient tending to zero as ft/ 3 ; note, however, that h a is not the 
distance which is usually encountered in the finite volume discretization of diffusion 
terms; consequently, the usual restriction for the mesh when diffusive terms are to 
be approximated by the two-points finite volume method (namelythe Delaunay 
condition) is not required here. We suppose that a > 1 and the convergence 
analysis uses < (3 < 2. 

Remark 4.1. At first glance, leaving the weight \p\ out, the stabilization term T s tat>,2 
may look as a Brezzi-Pitkaranta regularisation [2] , as used in [8] for stabilizing the 
colocated approximation of the Stokes problem, which would be rather puzzling 
since we use here an inf-sup stable pair of approximation spaces. However, using 
the equation of state (|1.2|) . we obtain: 

Tstab^-A 2 ( h K + h L f^ \pk+Pl\ (pk-Pl) 

a=K\L 

which shows, since A 2 = Ma 4 , that this term rapidly vanishes when approaching 
the incompressible limit. 

5. Existence of a solution and a priori estimates 

The existence of a solution to (|4.ip follows, with minor changes to cope with 
the diffusion stabilization term, from the theory developed in [10l section 2]. In 
this latter paper, it is obtained for fairly general equations of state by a topological 
degree argument. We only give here the obtained result, together with a proof of 
the a priori estimates verified by the solution, and we refer to [TU] for the proof of 
existence. 

Theorem 5.1. Let 6$ > and let T be a triangulation of the computational domain 
f2 such that 9 > 9q, where is defined by (|2.ip . Problem (|4.ip admits at least one 
solution (u,p) G Wh x L^; any possible solution satisfies px > 0, VK G T and: 

(5-1) ||u||i, 6 + \\p\\mn) + ||p||L»(n) + \p\r,p <C(f,M) 

where C(f,M) 6 K only depends on Q, A, f, M and 9q. 

Proof. Let (u,p) € Wh x Lt be a solution to (|4.1|) . Let pk — q{pk) for any 
K G T, and let p denote the vector (pk)keT- A natural ordering of the equations 
and unknowns in (|4.1b[) leads to a linear system of the form Mp = c, where c € R^, 
N is the number discretization cells, c £ M. N , c > 0, and where M is an M-matrix 
(in particular M^ 1 > and M~ l > 0). Therefore the i-th component of p reads 
Pi = M~ 1 c ■ ei — c • M~ t ei where ej is the i-th canonical basis vector of M. N . Since 
M ~ l > 0, we get M~ t ei > 0, and since M~ l ei 0, this proves that p. L > 0, which, 
in turns, yields pk > 0, Vif e T. Let us then prove the estimate (|5.ip . To this 
end, we take v = u in (|4. la|) and obtain: 

(5.2) / |Vhit| 2 da;— / pdiv^ttda;= / f udx. 

Jn Jn Jn 

Let us then multiply (|4.1b|) by A _1 [l + log(pic)] (see remark T5.2I below for an 
explanation of this choice) and we sum over K G T; dropping the terms which 
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vanish by conservativity, we then obtain: 
Ti + T 2 + T 3 = with: 

T l = A ~ X Y lo S(^) Y { V °,K PK - V a,K PL; 
KeT a=K\L 

T 2 = A- 1 h a Y \ K \ t 1 + 1 °S(PK)} [ P k - p*} , 
KeT 

T a = A- 1 ]T log(pK) Y ( h K + h L f ^{pk + Pl) {pk-Pl), 

KeT a=K\L a 

where the term \q(j>k) + q{pl)\ has been replaced by q(j>k) + q{jpl) m (|4.1b[) , thanks 
to the positivity of the pressure. Let us first write T± as: 

Tl = A^ 1 J2 Iog(px) J2 V <r,KP"> 
KeT a=K\L 

where p a is the upwind choice defined by (|4.2|) . Adding and substracting the same 
quantity, T\ equivalently reads: 

T 1 = A~ 1 Y;PK Y v^k + A- 1 ^ Y v a , K {po\og{p K ) - PK ). 

KeT a=K\L KeT a=K\L 

In the first summation, we recognize / p div^n da;. A reordering of the second 

Jn 

summation yields: 

Ti = / p div h u dx + A^ 1 Y] v^ K [{p„ log(p^) - p K ) - [p a log(p L ) - ph)\ ■ 
Jn =f . 

f Stint, 
a=K\L 

Pa = Pk = Pl if Pk = Pl, 
^p a log(pif) - p K = Pa log(pi) - pl otherwise. 
With this notation, we get: 

Tl = / p dlV h udx + A^ 1 V] y a,K {Pa - Pa) (log(p K ) ~ log(/5 L )). 

Jn ~T 

aetlnt, 
a=K\L 

In the last summation, we can, without loss of generality, choose the orientation 
of each edge in such a way that v a ,K > 0- With this convention, the term in the 
summation reads v 0i k {pk ~ Pa) (l°g(p/f ) — log(px,)), and is non-negative thanks to 
the fact that p a £ [min(px , p^), max(pK , Pl)] and the log function is increasing. 
We thus finally obtain: 



Let pa- be defined by 



Ti > p divftitda;. 



(5.3) 



Let us now turn to the estimate of T 2 . Since the function z >— > zlog(z) is convex 
for positive z and its derivative is z i— > 1 + log(z), we simply have: 



(5.4) 



T 2 > A- 1 h a Y \ R \ \P« lQ g(^) " P* lo g(^)] 



KeT 
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Finally, reordering the sums, the term T3 reads: 



T 3 = A- 1 (hK + h L f l -f (pk+pl) (j>k-Pl) (log(^) - log( Pi )) . 



a=K\L 



By concavity of the log function, we have: 



and thus: 



I logOif) - log(p L )| > 7 r \pk - Pl\ 

max(pK,PL) 



(5.5) T 3 > A~ x (h K + h L f 1 -^ (pk-Pl) 2 

a=K\L 



Summing equations (|5 . 2[) (j5 . 5[) and using Young's inequality, we obtain: 
||ti||i,6 +A- 1 / 2 \p\ T ,p <C(f,M). 

Furthermore, summing (|4.1bj) over K G T, we obtain that the mean value of the 
pressure p m is given by: 

p ™ = mL pdx=A ~ lp *- 

Using the inf-sup stability of the discretization, we get on the other hand: 

1 1 f 

\\P ~Pm\\h 2 (n) <— sup — j-j — r: / pdiv^da: 

Q vew h Fill, 6 Jn 

1 1 f 
= — sup I, / (Vfeit : V h v - f ■ v ) Ax. 

Q vew h lFl|i,6 Jn 

and the control of ||p||L 2 (fi) (or, equivalently, ||/o||L 2 (n) ) follows from the estimate 
for ||tt||i,6 ■ □ 

Remark 5.2 (On the choice of (log(px ))kgT as test function). At first glance, 
the choice of log(p#- ) to multiply (|4.1b|) in the preceding proof may seem rather 
puzzling. In fact, this computation is a particular case of the so-called "elastic 
potential identity", which is well-known in the continuous setting and is central 
in a priori estimates for the compressible Navier-Stokes equations [TBI HZ1 IS] • An 
analogous identity is proven at the discrete level, for the same discretization as 
here, in pjjl theorem 2.1]. 

For the particular case under consideration, an elementary explanation of this 
choice may be given. Indeed, it is crucial in the above proof that the quantity 
p a lies in the interval [vtm\{pK, Pl), niax(/9R-, pl)]- Let us suppose, without loss of 
generality, that < px < Pl and that, instead of the log function, the computation 
is performed with a non-specified increasing ond continuously differentiable function 
/; then we get: 

PL - Pk 
P ° f(PL)-f(p K y 
The condition p a > px is equivalent to: 

J_ > /(PL) - [(PK) 



Pk Pl - Pk 
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which is verified for / = log by concavity of the latter and, letting pl tend to pn , 
yields f'(x) <l/x. Conversely, the condition p a < pl yields: 

J_ < f(PL) - f(pi<) 
PL ~ PL - PK 

which, once again, is verified by the log function, and now implies f'(x) > 1/x. 

This limitation for the choice of the test function is the reason for the expression 
of the stabilizing diffusion term. 

6. Convergence analysis 
In this section, we prove the following convergence result. 

Theorem 6.1. Let (T^ m ^) m gN be a sequence of triangulations o/f2 such that h^ m ' 
tends to zero when m tends to +oo. Let us assume that this sequence is regular 
in the sense that there exists 9q > such that 9^ m ^ > 8q, Vm G N, where #( m ' is 
defined by (|2.ip . For m G N, we denote by Wfl and the discrete velocity 

and pressure spaces associated to 7"(" 1 ) and by (vS m \p( m ^) G W^™' ) x i[ the cor- 
responding solution to (|4.ip . with a > 1 and < [3 < 2. Then, up to a subsequence, 
the sequence (it^)meN strongly converges to a limit u in L 2 (f2) d and (p' m ') m gN 
converges to p weakly in ~L 2 (Q), where the pair (u,p) is a solution to in the 

following weak sense: 

u G Hj(fi) d , p G L 2 (n) and : 

(6.1a) f Vm : Vip dx - [ pdivipdx = [ f-ipdx, V^eC^Of, 

Jn Jn J q 

(6.1b) / pu-Vi/jdx = 0, V^GC^(fi), 

Jn 

(6.1c) f g(p) = M. 

Jn 

Proof. The proof is divided in three steps: we first show the existence of the limits 
u and p, then we pass to the limit in (|4.1aj) and (|4.1b|) . Since the equation of 
state is linear, the last equation is then a straightforward consequence of the weak 
convergence in L 2 (Q) of the (sub)sequence (p' m ') m£ M to p. 

Step 1: existence of a limit. 

By the a priori estimates of theorem 15.11 we know that: Vm G N, ||i^ m ^||i,6 < 
C(f,M). The compactness in L 2 (fl) d of the sequence (u( m )) me N, together with 
the fact that the limit u lies in Ho(S7) d , thus follows by applying theorem 13.31 to 
each component 1 < i < d. Once again by theorem 15.11 we have: Vm G 

N, ||p^||L 2 (n) < C{f,M). which is sufficient to ensure a weak convergence in 
L 2 (£!) of the sequence (p' m ') m£ N to p G L 2 (Q). 

Step 2: passing to the limit in (14. lap . 

Let ip be a function of C£° (0) d . We denote by i//" 1 -* the interpolation of ip in Wi" 1 ' , 
i.e. = r^tp where the operator r£ is defined by (|2.5p . Taking v — tp^ in 
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(|4.1a|) . we get: 

f V h U (m ^ : V h i/> (m) dx- f p {m) div h ip {m) dx = [ f- ip {m) dx, Vto e N. 
Jn Jn Jn 

Since the considered interpolation operator preserves the divergence (|2.6p . we have: 

/ p^ div/ l i/> (m " ) dx = / p^ m) divipdx — > / pdivipdx asm — > +00. 
Jn Jn Jn 



By the approximation properties of the interpolation operator (12. 8|) invoked com- 
ponent by component, we have: 



/ • ?// m - ) da; — > / f ■ ip dx asm — > 00. 

n 



Finally, we have: 

f V h tt (m) : V h ^ {m) dx = 

I V h u {m) : V h {ip [m) -ip)dx+ I V ft w (m ) : V^da; 
Jn Jn 



rp(m) rp{m) 
1 1 1 2 



Once again by (|2.8p , the term T[ m ^ obeys the following estimate: 

|T x (m) |< Wu^W^ H^-VHm <c(0 o )fc (ff,) |l« (m) lki. l</V(n) 



and thus tends to zero as m tends to +00. Integrating by parts over each control 
volume, the term T^™ 1 ^ reads: 

T 2 (m) = - / u (m) -Aipdx+Y^ I i u(m) ] V </> ' n ° d 7> 



where n CT is a normal vector to er, with the same orientation as that of the jump 
through tr. Applying Lemma 1 2. 41 for each component of Vi/>, a a being the relevant 
component of the normal vector n CT , we get: 



I [u (m) ] Vil> -n CT d 7 < c(9 )h^ \\u^\\ hb 1^1^ 

ate™ 



(Q) 



and thus tends to zero, while the first one tends to — J n u ■ Ai/> da; as m tends to 
+00. Since u G Hj(Sl) d , we may integrate by parts, and collecting the limits, we 
obtain (|6.1a|) . 

Step 3: passing to the limit in (|4.1b|) . 



Let ip be a function of C£°(il). Multiplying the second equation of (|4.ip by 



1/|X| / K ^(a;) da; and summing over K eT yields r 3 (m) +T 4 (m) +r^ m) = 0, Vm e N, 
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with: 



K<ET(">) 



rp{m) 



e ( e (p^+p^) (p^-pr)) 

• e r(™) \a=K\L " J 



where p^™ 1 ' is defined by (|4.2p and stands for the mean value of ip over K. 
Let q' m ) € Wj, be defined as q^ m \x) = c c-(™> ^i™ 3 p& ^(x), where ^ is 

the Crouzeix-Raviart basis function associated to a. The divergence of q^ is a 
piecewise constant function and reads: 

\/K e T(-), dW m > = j^j E <a P^ m) a.e. in K. 

' ' <r=A|L 



We thus have for T 3 



(m). 



T 3 (m) = E / ^W m3 dx 



Integrating by parts over each control volume, we get: 
T 3 (m) = - / W • <? (m) da; + E / ^ ' n " d7 

Vip-{p {m) u^)dx 



rp(m) 7 
J 6 



rp{m) 



By the respectively weak and strong convergence of (p' m ') m eN and («' m ')m£M to 
p and u in L 2 (f2) and L 2 (Sl) d , we have: 

/ Vtp • (p (m) w (m) )da; — > / V-0 • (pu) dx asm — > +oo. 
Jn Jn 

By the definition of q^ m \ the term Tg m ' reads: 
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with: 



r 8 (m) = £ ftp^l E «i7 ) ^(x)].n ff d 7 
= £ p<f) / ^ [« (m) ]-n CT d 7 

Since the integral of the jump of a Crouzeix-Raviart function over an internal edge 
of the mesh vanishes, the term Tg™^ can be estimated as follows: 



\Tt ] \<c^h^ £ pW / 



d 7 , 



where only depends on ip. Using the Cauchy-Schwarz inequality then yields: 

1/2 



\Tt l) \ <ci,hM £ M 1/2 pK/|[« (ro) ]| 2 d 7 ) 

<c^>(£ K\a\{p^f) l '\ £ ^/j[ uM ]| 2d ^ 

By the regularity of the mesh, the first summation is bounded by ||p^"'||l 2 (Q) 
while, by Lemma [2.21 the second one is bounded by c(0 o ) ||u (m) || 2 b . Let us now 
turn to the study of Tg . Since, for a' £ \ {er}, the integral of </v over a 
vanishes, and since the functions <f> a are bounded (namely \4>a\ < 1 in 2D, \<p a \ < 2 
in 3D) we get: 

/ V (^"P^MM*)] «i r r ) -n CT d 7 < C ^ ff H \p ( f-p^\ \u^\, 

J a 

where still only depends on ip. Since the function <p a i is non-zero over a = K\L 
only when a' belongs to the edges of K or L. only a limited number of terms are 
non-zero in Tg , in such a way that the difference — p ™ only involves two 
neigbouring cells or two cells sharing the same neighbour. Splitting the difference 
in this last case, using the previous inequality and the regularity of the mesh (in 
particular the fact that the ratio of the size of two neighbouring cells is bounded) 
and reordering the sums, we get for Tg™ 1 ' an estimate of the form: 

|T 9 (m) |<c £ hi\u^\ £ |/#°-pi m) |, 

(cr' = K\L) 

where the positive real number c only depends on ip and the regularity of the mesh 
and, thanks to this regularity, the set M a is such that a given edge K\L only appears 
in this sum a number of times bounded independently of m. Thus, thanks to the 
Cauchy-Schwarz inequality, we have: 

\Tt ] ?<c{ £ K\u^f) ( £ hi (pP-pP]* 

(a=K\L) 
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By the regularity of the mesh, the first term of this product is controlled by 
||ii( m )|| L 2(Q) and the second one by (ft/ m )) 2 ~ /3 \p( m ) \t,p ■ Consequently, thanks 

to estimate (|5.1|) . both Tg ' and Tg and thus also Tg tend to zero as m tends 
to +00, for any f3 < 2. 

Let us then examine the term X^™" 1 : 

T 7 m) = E / E (^ m) - pk } ) m*) • d* 

K£T(™) K <7=K\L 

Since V?/> is bounded in L°°(f2) , and since the functions <^> CT are bounded, we get: 

ir 7 (m) i <c„ e m E \pi m) -p { K n) \ l«4 m) l- 

Reordering the summations and using the Cauchy-Schwarz inequality yields: 
|T 7 (m) |<c„ E (\K\ + \L\)\p^-p[ m) \\u^\<^ (T^) 1 ' 2 (T^) 1 ' 2 , 

( CT =K|i) 

with: 

T& m) = E (|Ai + |i|)i4 m) i 2 



int 

(o-=if|L) 



Once again reordering the summation, we get: 



JfgT o-e£ K 



and thus, the term is controlled by \\u^ Hi^ro) , and T^^ is controlled by 

(/i( m )) 2 "^ |p ( - m - ) |r,/3 ■ By the a priori estimate (|5.ip . ^ thus tends to zero for any 
/3<2. 

We now turn to the terms T^™* and T^ m \ Since the sequence (p^)meN is bounded 
in L 2 (f2), the term J"!" 1 ' tends to zero for any a > 0. Reordering the summation in 

m(m) 

T 5 V ; , we get: 



rp(m) _ 
1 5 ~ 



(a=K\L) 



By regularity of ip, \ipK — V'il < c i/j (^a + ^l) an( l thus: 

lT^( m )l ^ II, 1 J. 17 ( ( m ) 1 ( m ) 

|T 5 I < c i> 2^ C 1 * + L) (Px + Pi 



(m) (m) 

Pk -Pi 



lot 

{<x=K\L) 
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Using the Cauchy-Schwarz inequality, we obtain: 
( 



|T 5 (m) |<c^/ 2 



1/2 



cre£. 



V ' — mt 
o=K\V) 



\p {m) W, (i 



J 



which, once again since the sequence (p^ m ^)meN is bounded in L 2 (f2), tends to zero 
by regularity of the mesh for (3 > 0. The proof is thus complete. □ 



7. Discussion 

To our knowledge, the convergence analysis performed in this paper seems to be 
the first result of this kind for the compressible Stokes problem (and, of course, more 
widely, for the compressible Navier-Stokes equations). Beside the convergence of the 
scheme, it also provides an existence result for solutions of the continuous problem, 
which could also be derived from the continuous existence theory ingredients for 
the steady Navier-Stokes equations (see [17] and references therein), but does not 
seem to be a direct consequence of the published literature: existence of strong 
solutions of the Navier-Stokes equations is known only for small data (e.g. [20] ) 
and existence of weak solutions is only proven for a particular class of equations of 
state (typically, p — p 1 with 7 > 3/2), this limitation being due to the presence of 
the convection term. 

A puzzling fact is that the present theory relies on two ingredients which are 
usually not present in actual implementations. Firstly, the stabilisation term T s tat>,2 
is needed in our proof to ensure the convergence of the discretization of the mass 
convection flux div(pu) and, to our knowledge, has never been introduced elsewhere. 
Secondly, the control of the pressure in L 2 (f2) relies on the stability of the discrete 
gradient (i.e. the satisfaction of the so-called discrete inf-sup condition), which is 
not verified by colocated discretizations; note that this argument is not needed for 
the stability of the scheme (see the proof of a priori estimates here and [TO] [7] for 
stability studies of schemes for the Navier-Stokes equations). Assessing the effective 
relevance of these requirements for the discretization should deserve more work in 
the future. 

An easy extension of this work consists in replacing the diffusive term —Ait in 
(|l.la|) by its complete expression — p Ait — ///3 V(dfvu) with p > (i.e. the usual 
form of the divergence of the shear stress tensor in a constant viscosity compress- 
ible flow) . Another less straightforward issue is the extension to more general state 
equations (for instance, p — p 1 with 7 > 1); it will be the topic of a further pa- 
per. Concerning higher order issues, let us note that the fact that the pressure 
is approximated by a piecewise constant function appears crucial in both stability 
and convergence proofs: extending this study to higher degree finite element dis- 
cretizations thus certainly represents a difficult task. Finally, let us remark that 
the present scheme relies on the approximation of the whole velocity vector at the 
interfaces. A less expensive scheme would be possible with a discretization it • n at 
the interfaces, as in the MAC scheme which is well known for the incompressible 
Navier-Stokes equations. However, such a discretization does not seem straightfor- 
ward on unstructured meshes. 
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